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It is believed that, much like a cat's cradle, the cytoskeleton can be thought of as a network of 
strings under tension. We show that both regular and random bond-disordered networks having 
bonds that buckle upon compression exhibit a variety of phase transitions as a function of tem- 
perature and extension. The results of self-consistent phonon calculations for the regular networks 
agree very well with computer simulations at finite temperature. The analytic theory also yields 
a rigidity onset (mechanical percolation) and the fraction of extended bonds for random networks. 
There is very good agreement with the simulations by Delaney et al. (Europhys. Lett. 2005). The 
mean field theory reveals a nontranslationally invariant phase with self-generated heterogeneity of 
tautness, representing "antiferroelasticity" . 



It is thought that the cytoskeleton consists of a poly- 
meric network which is maintained under tension [1, 2, 3], 
much like the child's game made of strings known as a 
"cat's cradle". Unlike a macroscopic "cat's cradle", the 
cytoskeleton is subject to fluctuating forces from both 
thermal motions and the noncquilibrium noise due to 
motor proteins and constant polymerization and depoly- 
merization events [4]. 

Upon freezing molecular fluids that resist compression 
undergo transitions to rigid solids. Here we show that, a 
caricature of the cytoskeleton, a model of interconnect- 
ing ropes that easily buckle also undergoes a variety of 
phase transitions. We study networks of ropes that bear 
tension when they are stretched beyond a critical thresh- 
old, but otherwise buckle. While the model's interactions 
contrast with the usual compression-bearing mechanical 
interaction between atoms or molecules, the "inverted" 
nonlincarity of ropes still leads to an onset of rigidity 
like that in crystals or glasses, but when the system is 
expanded rather than compressed. Depending on where 
we set the baseline of environmental pressure or tension, 
we may also say the system buckles upon compression. 

The cytoskeleton actually has both tension-bearing 
and compressional elements. The present theory, can be 
extended to include both elements. Elsewhere we have 
discussed the nonequilibrium aspects [5] that are also 
needed to capture the essential features of the tenseg- 
rity model [2, 3] of cell mechanics [6], but here we limit 
ourselves to the effects of equilibrium fluctuations. 

Delaney, Weaire, and Hutzler pioneered the study of 
such buckling lattices with neighboring sites that interact 
with an "inside-out" potential u(r; I) = 9(r— i)-|-(r— I] 2 , 
where 0() is the Heaviside unit step function [7]. The 
energy stored in a bond vanishes when the length of 
the bond r is less than a critical value I, but increases 
quadratically beyond the cutoff I. A two (three) dimen- 
sional network of bonds is thermo dynamically character- 
ized by its tension (negative pressure), as a function of 
its area (volume) — p2 = dF/ds (~P3 = dF/dv). 

In the low temperature limit, the entropy term can be 
ignored and elementary statics suffices to describe reg- 
ular networks. The square network, for example, has 



a tension as a function of the area expansion given by 
-p s 2 q (s/s ) = k(l - y/so/s) 9{s/s - 1). Rigidity only 
occurs for areas exceeding sq. Tensions for other lattices 
can also be computed p\ rl = 2\fip s 2 q and p 2 ex = p 2 "/3 
for the triangular and the hexagonal network respec- 
tively. For the 3D cu bic n etwork one finds at T = 0, 
-P s f(v/v ) = 2fc(l - ^7^) $/v^/v~ Q(v/v - 1) yielding 
a rigidity threshold at vq. We have set sq (or vq) as the 
area (or the volume) of the corresponding network con- 
structed from all bonds at the unit length. The tension 
of 3D networks has a maximum as a function of volume, 
because at large volume the energy density of the system 
decreases as u -1 / 3 . 

At high temperature or equivalently for soft springs, 
even without disorder, static models are insufficient and 
one always needs a statistical mechanical description. If 
the network is disordered, i.e., having a broad distribu- 
tion of onset length I, finding the tension is a nontriv- 
ial statistical problem even in the limit of low tempera- 
ture [7]. 

To analytically treat statistical cat's cradles, we adopt 
the self-consistent phonon technique, initially developed 
for anharmonic lattice solids with compression bearing 
interactions [8, 9]. In this approach one replaces the 
pairwisc interaction in the equilibrium density matrix 



cx P h/3££ ? 



j)] by a nearly equivalent effec- 



tive single-body term as IL, exp[— Hi(fi — Ri)]. Decoupling 
the manybody effects can be carried out by assuming a 
self-consistent one-body potential u has a quadratic form 
a{f-Rf. 

The self-consistent phonon method has not only been 
applied to crystals but also to noncrystalline solids [10], 
network glasses [11], and to problems of nonequilibrium 
structural assembly [12]. The self-consistent phonon 
method can be made the beginning of a systematic ex- 
pansion [8], but already has been shown to work very 
well for systems made up of compressive elements. 

We study both 2D and 3D systems. Long wavelength 
phonons lead to a diverging mean square displacement 
(MSD) with increasing system size for 2D systems, in a 
strict sense. However this divergence is weak (logarith- 
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mic) and there is still a crossover reflecting rather rapid 
decrease in fluctuations of the mean square nearest neigh- 
bor distance upon extension. In fact, the self-consistent 
phonon method, when used for the 2D system predicts 
the correct onset of rigidity for the random 2D network 
at T = as we shall see. 

In the self-consistent phonon theory the partition func- 
tion z(£) = [J dr exp(— e(£) — ct(£,)r 2 )] N is evaluated as a 
function of expansion ratio £ = s/sq or v/vq. Both the 
on-site free energy e and the phonon frequency a depend 
on the ratio of the mean neighbor distance R to the crit- 
ical extension /, R/l = We can use the chain rule 

to calculate tension, e.g., -p 2 = §7 = -^0^-- 

In three dimensions the Mayer /-bond v(R; I, a) for 
the interaction u(r; /) when averaged over the Gaussian 
fluctuation of a neighboring particle around its mean po- 
sition R having its own spring constant a gives 

exp[-pv(R)] = (^)-S J r 2 dr J dfi e -a(i'-j*) a -0«M 

This can be expressed using the error function erf(). The 
2D averaged Mayer /-bond can be expressed as a numer- 
ical integral of modified Bessel function Iq. 

An interesting property that can be calculated within 
the self-consistent phonon theory is the fraction of the 
bonds bearing force (that are taut), i.e., having a length 
larger than the intrinsic cutoff I. These are functions of 
a and R: qi = 2ae~ aR rdrIo(2aRr)e~ ar and 93 = 

l + ^[(R-l)^-^[(R+l)^+ ia % 1/2 [e^ R -^ 2 - 
e~ a ( R + 1 ) 2 } for 2D and 3D cases. 

The self-consistent phonon frequency a for several dif- 
ferent regular networks is shown in the left panels of 
Fig. 1 . Unlike the situation for the systems made of hard- 
spheres [5], there is always some nonzero self-consistent 
solution a. Entropy always leads to some rigidity. There 
is, however, a rapid crossover from low a to high a ac- 
companying the expansion of the system. Generally ei- 
ther high coordination number z or larger microscopic 
rigidity-temperature ratio b (= k/3/2) lead to a larger a. 
The fraction of "loose" bonds, (1 — q) predicted from the 
phonon theory is plotted in the right panels of Fig. 1. 
This fraction decreases as the network expands. 

We have shown there is a distinct bifurcation of a 
that also occurs when the volume ratio £ becomes small 
enough for 3D networks with high rigidity-temperature 
ratio b, This bifurcation appears in all three types of 3D 
networks studied here. The fee networks most clearly 
display this bifurcation, followed by the bec, and then 
the sc networks at smaller £ and higher b. As shown 
in Fig. 2, usually there is only one stable solution of the 
self-consistent equation a ta gged = f (^neighbor) ■ However, 
under certain conditions (small volume and large rigidity- 
temperature ratio, or pictorially, when the slope at the 
intercept is not in the range of between -1 and 1), the 
self-consistent solution becomes unstable, and another 
pair of solutions a,\ and ah satisfying ai = fioth) and 
a,h = f(ai) emerges as the stable solution. This stable 



pair solution indicates the possibility of having a non- 
translationally invariant solution, in which a low a par- 
ticle is surrounded by high a neighbors and vice versa. 
The resulting spatial heterogeneity is self-generated. In 
the bec and sc networks, this solution can be easily in- 
terpreted physically: there is a symmetry-broken phase 
of bec or sc networks that can be thought as an "antifer- 
romagnctic" (AF) one with the underlying lattice being 
split to two subnetworks (sc — > fee or bec — > sc): one of 
the two equal subnetworks bears more tension than the 
other. We call this the "antiferroelastic" phase. In the 
fee case, since there are four particles in the unit cell, one 
high a particle cannot be completely surrounded by low 
a particles. Nevertheless one can have two possible or- 
dered AF ground states similar to the crystal structures 
of AuCu and AuCu3. To test this possibility, we set up 
a more elaborate self-consistent solution having four val- 
ues of a = (a ,a x ,a y ,a z ) at (000), (Oil), (101), (110) 
respectively, a = f(a x ,a y ,a z ), a x = f(a ,a z ,a y ), 
at y = f{a z ,a ,a x ), a z = f(a y ,a x ,a ). By explicitly la- 
beling all four particles in a unit cell (as done similarly 
in the AF Ising fee system [13]), we searched but failed 
to find these explicit symmetry-broken solutions. Thus 
we believe this frustration of fee network would lead to 
a random phase which is analogous to a self-generated 
glass. This possibility may potentially be resolved using 
a replica calculation [14, 15]. 

We checked our self-consistent phonon theory using ex- 
plicit computer simulations of the network with thermal 
noise. We show the results for the bec case in Fig. 3. The 
simulations were performed with stochastic dynamics on 
n 3 unit cells with periodic boundary condition for vari- 
ous rigidity-temperature ratios b and nonlincarity param- 
eters / for various values of n. We see the self-consistent 
phonon theory agrees quite well with the simulation. For 
the simulation of size 5 , the errors for the MSD of the 
particle position between theory and simulation are only 
about 1%. We extrapolated to the thermodynamic limit 
based on a series of simulation of different sizes of net- 
works. A Similar good agreement is achieved for 2D net- 
works when the MSD of relative neighbor position from 
simulation is compared with theory. 

To check the prediction of spontaneous symmetry 
breaking of neighboring a values requires a large ratio 
of l/R ~ 10 necessitating simulation using a very large 
system at very low temperature. We estimate the neces- 
sary system size to be at least 10 6 . We were unable to 
carry out such a large scale simulation. 

We plot the tension vs. extension in Fig. 4 (a) and 
(b). Above a critical tension, there are two solutions 
for p(s/so) = P- The solution with dp/ds < is ther- 
modynamically stable for a constant pressure ensemble. 
The system always appears to be stable in two dimen- 
sions. For 3D networks, As shown in Fig 4(b), there is 
single (unstable) solution at high temperature for con- 
stant P3; but there are three (one stable and two un- 
stable) solutions at low temperature. The values of the 
critical rigidity-temperature ratio for this transition are 
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FIG. 1: The phonon frequencies a are shown in the left 
panels as functions of the area (a) or volume (b) of networks 
for various b = f3k/2, and as functions of the linear expansion 
parameter with constant b — 1 for various types of networks 
(c). The fractions of free bonds 1 — q of corresponding setups 
are shown in the right panels (d)-(f). Note that (a), (b), (d), 
and (e) share a legend while (c) and (f ) share the other legend. 




FIG. 2: The illustration of the onset of bifurcation solutions 
for a specific case. 



b* = 1.04, 0.76, and 0.48 for sc, bcc, and fee networks 
respectively. Thus the network is intrinsically unstable 
for the constant p ensemble for high temperature or low 
spring constant, i.e. b < b* , regardless of supplied ten- 
sion —p. For b > b*, there is a critical tension — p* above 
which the system is no longer stable. Generally a larger 
b implies a larger stable region. 

In the large s/sq or v/vq limit, the Dulong-Pctit law 
provides pressure that can be based on dimension count- 
ing alone, giving p 2 — p° = NkBT/{2S) and pz — p° = 
Nk B T/(2V). Here p° = p(T = 0) is the pressure from 
the mechanical static calculations at zero temperature. 
The factor | arises because only the potential energy is 
counted here, the kinetic energy counts the other half. 
Thus for large expansion ratios, pi — > const, and p^, — > 0, 
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FIG. 3: The comparison of the MSD of bcc networks between 
theory and simulation for various parameters R = y/3/2 x A, 
b, and I. The corresponding results for square networks are 
shown in the inset (R = A). 



which are consistent with the tension curves shown in 
Fig. 4(a) and (b). 

So far we have focussed on the effects of thermal dis- 
order alone. We now discuss the effects of bond disor- 
der. For illustration, we concentrate on the 2D random 
network, especially in the (3k — > oo limit. Following De- 
laney et al. [7], we studied a "flat" random distribution 
of string lengths. More general distributions present no 
new challenges to our method. 

As shown in Fig. 4(c) for the triangular network, 
phonon frequencies a rise quickly with expansion at zero 
temperature when b = oo. The smallest area s* /sq sat- 
isfying a _1 (s/sn) = is the critical value for the onset 
of macroscopic rigidity (mechanical percolation). Since 
the spring constant k and temperature are always bun- 
dled together in the treatment, once the system is taut, 
a — > oo for 1/6 = 0. For finite temperature, a is always 
finite. Th rigidity onset s* /sq calculated here quanti- 
tatively agrees with the simulation results in Ref. [7] as 
a function of width of distribution of string length w. 
As expected, a larger w means a lower a generally. We 
also display in Fig. 4 the fraction of vertices with given 
number of taut bonds P n (d) and the mean fraction of 
taut bonds q (c). These properties have also been ob- 
tained by the simulations of Delaney et al. We can only 
make a qualitative comparison for these two properties 
since our calculations were performed with finite b, cor- 
responding to T 7^ while the simulations were strictly 
done at T — 0. Our theory does not allow us to calculate 
the loose bond fraction at T = due to the bundling of (3 
and k. An approximation was made for the calculation of 
p(n), the fraction of vertices with given number n bonds 
taut based on the mean fraction of taut bonds q, i.e., 
P( n ) = n \(z-ny. Q n 0- — q) Z n - Here Z is the coordination 
number. 

As shown in Fig. 4(f), for bond-disordered 3D net- 
works, again there is a rigidity onset, which is similar to 
those of 2D networks. Among the cases studied only the 
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FIG. 4: The applied tensions as functions of area expansion of 2D triangular network for various b are shown in (a) with 
the zoom-in at the small s/so region shown as an inset. We show the tensions of 3D networks in (b). The random bond 2D 
triangular network results are shown in (c)-(e). The effects of the random bond distribution w on the phonon frequency a are 
shown in (c) with the inversion 1/a shown in the inset of (c) to better illustrate the critical area for onset of rigidity. The 
portion of vertices with n bonds tauted P„ and the mean taut ratio q are shown in (d) and (e) respectively. We show the 
inversion 1/a as functions of volume for various 3D random networks in (f). 



w = 1 fee network showed the antiferroelastic structural roelasticity. This work was supported in part by NSF- 
inhomogeneity bifurcation. The self-consistent phonon CTBP. 
method suggests random bond-disorder destroys antifer- 
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